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MODELING  AND  INFERENCE  FOR  POSITIVELY  DEPENDENT 


VARIABLES  IN  DICHOTOMOUS  EXPERIMENTS 


ABSTRACT 

Multivariate  models  with  positively  correlated  components  have  found 
wide  applicability  in  reliability  and  biostatistics.  Perhaps  the  best 
known  and  most  widely  used  such  model  is  the  multivariate  exponential 
distribution  due  to  Marshall  and  Olkin  (JASA,  1967) .  We  study  a  discrete 
analogue  of  the  latter  model.  Specifically,  we  consider  a  model  for 
random  vectors  Y  whose  components  are  positively  correlated  and  have 
Bernoulli  marginal  distributions.  The  construction  of  the  model 
reflects  the  fact  that  the  k  component  system  under  study  may  be 
subjected  to  independent  shocks  selectively  fatal  to  any  subset  of 
components.  A  special  representation  of  the  probability  function  of 
Y^  is  developed  which  proves  useful  in  the  inference  questions  pursued. 
While  maximum  likelihood  estimation  of  the  model  parameters  proves 
intractable,  we  obtain  in  closed  form  an  alternative  estimator  which  we 
show  to  be  asymptotically  equivalent  to  the  MLE  and,  in  fact,  equals 
the  MLE  with  limiting  probability  one.  Similar  results  are  obtained 
for  a  natural  submodel  whose  parameter  space  is  of  substantially  lower 
dimension.  A  Monte  Carlo  study  sheds  light  on  the  sample  size  needed  . 
for  the  asymptotic  results  to  take  hold.  j  / 
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I.  INTRODUCTION 
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Consider  the  vectors  in  f 0 , 1}  as  being  ordered  from  the  smallest, 
(0,0,0, ... ,0)  to  the  largest  (1, 1,1, . . . ,1)  according  to  the  usual  order¬ 
ing,  that  is,  lexicographically.  Attach  to  the  itl1  such  vector  the 
probability  p^  €[0,l],  for  1*1,2, ... ,2^,  subject  to  the  constraint 

2k 

J  p,  ■  1.  The  description  above  defines  the  most  general  multivariate 
i-1  1 

Bernoulli  distribution,  and  any  k-variate  distribution  with  Bernoulli 
marginals  can  be  fully  described  by  identifying  the  vector  p  above. 

A  useful  reference  on  such  distributions  is  the  paper  by  Bahadur  (1961). 
Inference  questions  for  the  general  multivariate  Bernoulli  distribution 
are  easily  resolved;  for  example,  the  maximum  likelihood  estimate  of  pa 
from  a  sample  of  size  n  is  simply  the  relative  frequency  of  occurence 
of  the  i  vector  in  {0,1}  .  There  are,  however,  a  number  of  interesting 
subclasses  of  multivariate  Bernoulli  distributions  which  arise  naturally 
in  applications  but  lend  themselves  less  readily  to  statistical  analysis. 
This  paper  is  dedicated  to  the  study  of  one  such  class. 

Multivariate  models  which  postulate  positive  dependence  among  the 
components  of  a  random  vector  have  found  substantial  applicability  in 
reliability  theory  and  in  biostatistics.  For  example,  the  lifetimes  of 
items  on  test  may  well  be  positively  correlated,  particularly  when  these 
items  are  subject  to  shocks  that  threaten  two  or  more  items  simultaneously, 
or  when  some  items  experience  wear  which  serves  to  increase  the  load  on 
other  items  on  test.  The  multivariate  exponential  distribution  of 
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Marshall  and  Olkln  (1967)  arises  naturally  in  several  reliability  con¬ 
texts,  Including  a  scenario  Involving  a  collection  of  shocks  selectively 
fatal  to  one  or  more  components  In  a  coherent  system.  In  constructing 
competing  risk  models  In  problems  arising  In  the  health  sciences,  one 
often  encounters  positive  dependence  due  to  the  reduction  of  resistance 
to  one  disease  that  might  be  caused  by  the  presence  of  a  second  disease. 
The  model  to  be  studied  here  may  be  viewed  as  a  discrete  analogue  to  the 
multivariate  exponential  distribution  in  that  we  will  motivate  it  in 
much  the  same  way.  As  we  shall  see,  however,  the  multivariate  Bernoulli 
model  we  study  is  quite  general  and  may  be  viewed  as  a  discretized 
version  of  any  continuous  multivariate  model  describing  the  joint  life¬ 
time  of  several  components  subjected  to  selectively  fatal  shocks.  Al¬ 
ternative  notions  of  positive  dependence  in  reliability  have  been 
studied  by  several  authors,  including  Esary,  et  al.  (1972)  and  Shaked 
(1975).  It  is  well  known  that  these  notions  arc  all  equivalent  when 
dealing  with  Bernoulli  variables 

Consider  a  k-component  system  whose  lifetime  is  under  study.  Sup¬ 
pose  the  status  of  the  system  is  to  be  observed  at  some  distinguished 
time  Tq  which  could,  for  example,  be  the  so-called  mission  time  of  the 
aystem  or  of  the  individual  components.  The  observation  may  be  recorded 
as  a  vector  Y  ■  (Y^Yj, . . .  ,Yk)  where,  for  each  i  ,  Y^  la  one  or  aero 
depending  upon  whether  the  1C^  component  is  working  or  has  failed  by 
tine  Tq.  We  will  formulate  a  model  for  the  distribution  of  Y  which 
reflects  the  fact  that  the  system  may  be  subjected  to  independent  shocks 


3 


selectively  fatal  to  any  subset  of  the  k  components.  To  this  end, 

define  “  (0,l)  -  (0),  that  is,  let  be  the  collection  of  all  k-tuples 

of  zeros  and  ones  with  the  exception  of  the  zero  vector.  Each  element 

s€*^  corresponds  to  a  shock  selectively  fatal  to  those  components  i  for 

which  s.  -  1.  For  each  s€«/,  ,  define 
i  k 


Z 

s 


if  the  shock  corresponding 
to  s  occurs  by  time  Tq. 

otherwise. 


(1.1) 


and  define 


P 


8 


1). 


(1.2) 


(Our  notation  is  precisely  that  of  Marshall  and  Olkin  (1967).)  Denoting 
the  binomial  distribution  with  parameters  n  and  p  by  3 (n,p) ,  we  have 
that  Z  ~B(l,p  )  for  each  s€*/  .  We  may  now  define  the  aforementioned 

8  8  '*•  K 

vector  Y  componentwise  as 


Y 


i 


(b€ 


i 


min  Z  . 

(s€  J.  :  s .“l)  ~ 

~  k  i 


(1.3) 


If 

We  will  use  the  notation  MVB(2  -1)  to  denote  the  distribution  of  Y  in 
equation  (1.3),  and,  In  general,  the  notation  MVB(n),  for  various  Integers 
n  ,  to  denote  particular  submodels  of  MVB(2  -1)  with  exactly  n  para* 
meters.  Equation  (1.3)  may  be  Interpreted  to  mean  that  the  i^  component 
will  survive  until  time  Tq  If  and  only  if  no  shock  that  is  fatal  to  the 
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it*'  component  occurs  by  time  Tg.  The  shock  that  is  fatal  to  the  ith 
component  alone  has  the  dual  interpretation  of  being  the  embodiment  of 
all  factors  that  might  prove  fatal  to  the  i^  component  but  have  no 
effect  whatever  on  any  other  component. 

In  section  II,  we  study  the  properties  of  the  shock  model  we  have 
described  in  the  preceding  paragraph.  In  particular,  a  representation 
of  the  probability  mass  function  of  Y  (or  any  subvector  of  Y  )  is 
developed.  This  representation  proves  to  be  quite  useful  in  the  statis¬ 
tical  Inference  we  pursue  in  later  sections.  In  section  III,  we  consider 
maximum  likelihood  estimation  of  the  parameters  of  our  multivariate 
Bernoulli  distribution.  We  Investigate  and  comment  on  the  difficulties 
involved  in  obtaining  the  MLE  in  closed  form,  and  we  then  produce  (in 
closed  form)  an  alternative  estimator  which  we  show  is  asymptotically 
equivalent  to  the  MLE.  In  fact,  the  proposed  estimator  is  equal  to  the 

MLE  with  limiting  probability  one.  In  section  IV,  we  investigate  maxi- 

fc. 

mum  likelihood  estimation  for  a  natural  submodel  of  MVB(2  -1) , 
and  derive  results  which  are  comparable  to  those  in  section  III.  In 
section  V,  we  summarize  the  results  of  a  modest  Monte  Carlo  study  which 
sheds  light  on  the  sample  size  required  to  get  reasonable  results.  The 
simulation  study  motivates  our  discussion  of  the  primary  domain  of 
applicability  of  our  results,  namely,  that  of  systems  consisting  of 
components  subject  to  "rare  shocks."  We  summarize  our  results,  indicate 
some  promising  directions  for  future  work,  and  make  a  number  of  con¬ 
cluding  remarks  in  section  VI. 
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II.  CHARACTERISTICS  OF  THE  MVB(2k-l)  MODEL 

To  motivate  the  distribution  theory  presented  in  the  Lemma  and 

Theorem  below,  we  begin  by  examining  the  simplest  model,  that  for  a 

two-component  system.  When  k*2,  we  are  concerned  with  the  representation 

2 

of  P(Y  ■  y)  for  y  €(0,l)  in  terms  of  the  three  model  parameters 

P01’  P10  *n<*  pu*  Clearly,  the  following  representation  is  one  possibility: 

pq  -  (i.i))  -  P10P0iPn 
pot  -  (o,d)  -  a-p10>p0iPii 

(2.1) 

P(Y  -  (1,0))  -  (l-P0l)P10Pn 

P(Y  -  (0,0))  -  (l-pu)  +  Pll(l-p10)(l-p01). 

The  character  of  such  representations  for  arbitrary  k  Is  already  apparent 

from  equation  (2.1).  ■  It  would  seem  that  in  general  the  probability 

P(Y  ■  y)  could  be  written  as  a  sum  of  terms,  each  term  being  a  product 

of  certain  parameters  p  or  their  complements  (1-p  ).  A  systematic 
— "  8  8 

approach  to  this  sort  of  representation  would  proceed  as  follows:  Con¬ 
sider  the  tree  consisting  of  the  eight  branches  which  represent  all 
possible  combinations  of  the  three  dichotomies  (Z  “0  or  Z  -l[s€*^,,}. 

This  tree  Is  pictured  below: 
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Probability  Y 


P01P10P11 

(1,1) 

P01P10<1“P11) 

(0,0) 

P01(1-P10^P11 

(1,0) 

p01(1'p10)<1-pll) 

(0,0) 

(1-p01)p10pll 

(0,1) 

<1“p01)p10(1-pll) 

(0,0) 

(l-p01)(l-p10)pu 

(0,0) 

(1-p01)(1"p10)(1"pll) 

(0,0) 

(2.2) 


The  probabilities  associated  with  each  of  the  four  possible  values  of 
Y  may  be  obtained  from  this  tree  simply  by  grouping  branch  probabilities 
appropriately.  A  general  representation  for  the  probability  mass  function 

ft 

of  Y  ~  MVB(2  -1)  may  be  obtained  In  the  same  manner.  However,  since 

It 

there  are  2  -1  shocks,  the  tree  from  which  probabilities  are  to  be 

2^-1 

Identified  will  have  2  branches.  Thus,  the  representation  along 
these  lines,  while  conceptually  simple,  poses  some  practical  difficulties. 
This  approach  can  be  made  more  appealing  by  applying  combinatorial  ar¬ 
guments  which  facilitate  the  counting  of  ways  in  which  a  specific  value 
of  Y  can  arise.  We  do  not  pursue  thia  further,  however,  because  we  have 
found  this  type  of  representation  less  useful  in  inference  problems 
than  the  representation  described  in  Theorem  2.1.  As  a  closing  remark 
on  the  representation  discussed  above,  we  note  chat  the  probabilities 
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2-1 

associated  with  the  2  branches  of  the  aforementioned  tree  are  simply 
the  terms  of  the  polynomial  expansion  of  the  right  hand  side  of  the  equation 
below 

1  -  TT  (P  +  <1-P  ))  (2.3) 

{b6^}  *  ~ 

-  k 

Before  establishing  an  alternative  representation  of  P(Y  =  y)  in 
the  general  model,  we  introduce  the  following  notation,  where  all  vectors 
are  k-dimensional. 


0  =  (0,0,0, ... ,0) 

J!  =  U  (0)  (2.4) 

k  k 

1  =  (1,1,1, ... ,1) 

y*  =  1  -  y  (2.5) 

~  ~  k 

I y !  =  | <y,*  —  > I  =  I  |yj  (2.6) 

1  *  i-i  1 

xy  -  (x1y1,x2y2,...,xkyk)  (2.7) 

!k(y)  -  (s  yt“l  *  s^l,  i-1 . k)  (2.8) 


Now,  let  j  be  an  integer  less  than  or  equal  to  k  ,  and  let  N^-fn^  i-l,...,j) 
be  a  set  of  j  integers  such  that  1  £  n^  <  n2  <• ■ .<  n^  £  k.  We  will  have 
occasion  to  consider  mappings  nN  :  +  of  the  form 


v 5  yvv- 


). 


(2.9) 
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When  confusion  seems  unlikely,  the  mapping  TT 

J 

or  simply  by  TT, 

We  need  the  following  preliminary  result. 


will  be  denoted  by  TT 


J 


Lemma  2.1.  For  any  fixed  but  arbitrary  set  of  j  ordered  positive 
Integers,  with  j  <  k  and  n^  <  k. 


P(Yn  "Yn  -1> 

n,  n„  n. 


■  TT  p„ , 


€  A 


(2.10) 


where  A  -  [s  €/,  :s  "l  for  some  i-l,...,j} 
~  k  n  > 


Proof.  P(Y  -  . . .  -  Y  -1) 
nl  "j 


P(  TTz  Tt  Z  “  l) 

8  *1  *  S  *1  * 


‘j 


-  P<z8  ’  i 


9  s  *1  for  some  i  ) 
ni 


-  TTp(zb  -  1). 

5€a  A 


the  last  step  being  a  consequence  of  the  independence  of  the  Z  variables. 

Theorem  2.1.  Let  be  a  fixed  but  arbitrary  set  of  j  ordered  positive 
integers,  with  J  <  k  and  n^  <  k.  For  every  we  have 

P(TTY-j)  -  S  <-l)'~'^  TTp.  (2.1D 

*€lj(*>  Nr,,,- 

where  tt  ■  ttn  «nd  An>w  -  :  (TTs)  w  *  0}. 


j 
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Proof :  The  proof  proceeds  by  Induction  on  |y*l,  the  number  of  zero  components 

of  the  vector  y  ,  where  the  dimension  J  of  y  may  vary  between  1  and  k  . 

We  show  that  for  each  fixed  k,  the  representation  in  (2.11)  holds  for 
each  J«l,...,k  and  for  all  values  of  jy*|  <  k.  If  |y*|  -  0,  (2.11) 
follows  from  Lemma  2.1.  Now  suppose  that  |y*|  ■  h,  where  1  <  11  5 
Assume  that  for  each  J“l,...,k,  and  for  all  possible  mappings  tt^  of  the 
form  (2.9),  the  representation  in  (2.11)  holds  for  P(n^  y  *  £ )  whenever 
s€*^  is  such  that  Js*l  <  h.  Since  y  has  at  least  one  zero,  we  may 

~  J  ~ 

*  t  » 

assume  without  loss  of  generality  that  y^  »  0.  Define  TT  •  J ^  -*  J  ^  by 
★  * 

TT  s  *  TT  "  (92  9  S3  *  *  *  #  *  8  j  ^  9 

and  define  by 

~  J 


if  i  >  1 


if  l  -  1. 


We  then  have 


P(TT*TT  Y  -  TT*y)  -  P(Y  -y . Y  -  y .) 

J  ~  ~  n2  1  nj  J 

-  P(TTY  -  y)  +  P(TTY  «  y°) 


P(ttY-y)  -  P(n*!TY-TT*y)  -  P(TTY-y°). 


Since  both  TT*y  and  y°  have  exactly  h-l  zero  components,  we  have  by 


the  induction  hypothesis  that 
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The  representation  in  (2.11)  gives  the  distribution  of  the  random 
vector  Y  and  shows  that  the  marginal  distribution  of  any  subvector  has 
the  same  form.  We  state  as  a  corollary  the  representation  result  to  be 
used  in  our  study  of  inference  questions. 

Corollary  2.1.  Let  Y  MVBCZ^-l)  for  k  >  2.  Then  for  any  y  g  J'  , 

~~  k 

P<2-£)  -  S  (-l)W  TT  P8  (2.13) 

w€ik(,x)  Ue*V  ~ 

8  W  4  0  } 

A/  J 


We  have  mentioned  above  that  the  model  under  discussion  has  the  property 
that  components  of  a  vector  Y  described  by  the  model  are  positively  corre¬ 
lated.  This  is  easily  demonstrated  by  noting  that 


Y^  ~  B(l,  TT  p^),  so  that 

r*i, J  k  r 

C"<VV*  TT  •(  TT  ».X  TT  O^0* 

r.V  Jt~6’,k!vl5 


with  equality  holding  if  and  only  if  pg  •  0  for  acme 

or  p  •  1  for  all  a  for  which  s.  »  s.  •  1. 
a  ~  *  1  J 


s  €  u  {s  ;  •  »l) 

~  r«i,J~  k  r 
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The  simplest  possible  multivariate  Bernoulli  distribution  is  the 
model  in  which  components  Y^  are  independent  and  identically  distributed 
as  B(l,p).  Conditions  on  the  parameters  of  the  model  MVB(2k-l)  which 
result  in  partial  or  complete  simulation  of  the  model  with  i.i.d.  com¬ 
ponents  are  as  follows:  (i)  Dependent  components  with  identical  marginal 

distributions  result  if  and  only  if  p  ■  p  ,  for  any  s,  s'  for  which 

s  s  ~  ~  k 

\s\  *  |s,| ,  (ii)  Independent  components  with  nonidentical,  nondegenerate  * 
marginal  distributions  result  if  and  only  if  p  ■  1  for  every  sf/ 

S  ~  K 

for  which  |s|  >  1  (that  is,  only  shocks  fatal  to  a  single  component  have 
a  positive  probability  of  occurrence),  (iii)  The  Y^,  i*l . k  are  non¬ 

degenerate,  independent  and  Identically  distributed  iff  the  conditions 
on  s  €/.  in  (1)  and  (ii)  above  are  simultaneously  satisfied. 


III.  ESTIMATION  FOR  MVB(2  -1). 


We  begin  this  section  with  an  examination  of  the  problem  of  maximum 

likelihood  estimation  of  the  parameter  for  the  model  MVB(2  -1).  To 

fix  ideas,  we  briefly  discuss  the  case  k-2,  that  is,  the  two  component 

2 

model.  Let  Y,,...,Y  be  a  sample  of  size  n  from  MVB(2  -1),  and  for 

each  y  €  [o, l]  ,  let  N  be  the  frequency  of  the  observation  Y  *  y  .  A 
^  y  m  ^ 

complete  discussion  of  the  case  k«2  Involves  consideration  of  the  15  data 
configurations  where  each  N  is  zero  or  nonzero,  with  at  least  one  being 
nonzero.  Suppose,  for  example,  N  >  0  for  each  y  .  Then  the  likelihood 

y  ~ 

function  may  be  written  as 


L  -  [(1 -?„)+?„<! -plft)(l-pn1)3  °°[(1-p0)p10pu] 


10 


'11'  ■  Kll’ 


10/v*  l'01" 
N. 


N, 


•Co.-  Pio>p01pll-^  ^P10P01PU^ 


(3.1) 


13 


It  is  clear  even  in  the  simple  case  being  considered  that  maximum  likeli¬ 
hood  estimation  is  fairly  cumbersome.  The  system  of  equations  to  be 
solved,  that  is,  L  *  0  for  a€^2,  are  highly  nonlinear  and  each 

involves  all  parameters  in  a  complex  way.  it  can  be  shown  (by  direct 
maximization)  that  the  MLE  when  N  >  0  Vy  is  given  by 


:  .  "li 

10  N11  +  N01 


_  -  11 
01  N11  +  N10 


(3.2) 


when  NnNoo  *  NoiNio’  and  t8 


n  p10p01 


Nn  +  Nio 


_  Nll  +  N01 


(3.3) 


P11  “  1 


when  NuNQ0  <  N01N1().  We  are,  in  fact,  able  to  find  the  MLE  in  closed 
form  for  the  case  k-2,  but  will  not  pursue  direct  maximum  likelihood 
estimation  further  for  several  reasons.  It  is  clear  that  direct  maximum 
likelihood  estimation  is^unpromising,  involving,  in  the  case  of  arbitrary 
k,  the  examination  of  22  -l  data  configurations,  each  associated  with  a 
highly  complex  system  of  (2k-l)  equations.  Moreover,  there  is  no  guarantee 
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that  any  one  of  these  systems  has  a  closed  form  solution.  Indeed,  we 
have  found  that  there  is  no  closed  form  solution  of  the  likelihood 
equations  in  some  of  the  special  cases  we  have  studied.  We  will 
thus  turn  our  attention  to  the  construction  of  an  alternative  estimator. 
We  show  in  the  sequel  that  this  estimator  is  asymptotically  optimal  and, 
in  fact,  is  equal  to  the  MLE  with  limiting  probability  one. 

Jq 

Let  Y. . Y  be  a  randcm  sample  from  MVB(2  *1)»  and  define  N 

/w  »vn  y 

for  each  y€»V  as 
k 


frequency  of  occurrence  of  Y  =  y  . 


If  Q  s  P(Y*y),  then  the  vector  N  *  (N..,...,!*-),  with  components 
y  u  l 

rw  ^  a# 

lexicographically  ordered,  has  a  multinomial  distribution  with  parameters 
n  and  Q  .  The  MLE  of  Q  is  given  by 


for  l €  ''k  * 


(3.4) 


The  estimator  we  develop  is  derived  from  an  attempt  to  invert  the 
functional  relationship 


Q  - 

y 


2 

w€  I,(y) 


1  w  -  y  1 
(-1)  ~  ~ 


TT 


{s  :  s  w  i  0} 

9SS  K  ~ 


(3.5) 


for  each  y€»^  to  find  p  as  a  function  of  Q.  Under  well-known 
conditions,  the  estimate  of  p  obtained  from  (3.5)  and  (3.4)  will  be 
the  MLE.  These  conditions  do  not  obtain  in  the  present  problem;  we  are 
nevertheless  able  to  use  the  invariance  property  of  MLE's  to  advantage. 
We  proceed  with  the  inversion  of  the  relationship  in  (3.5). 
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Let  yC*'.  be  such  that  jy|  “  1.  Then  |y*|  -  k-1,  and  we  obtain 


from  (3.5) 


Q  *  -  TT  Ps  -  TT  P 

Z  a  t  y  ~  sf/ 


s 


(3.6) 


Since  ■  TT  pg,  we  obtain  from  (3.6)  that,  provided  #  o, 


~  i^k 


Q  * 

P  -r jl  +ii 

Z  L  Qi  J 


Q,  +  Q  * 


(3.7) 


Now  let  y€Vv  be  such  that  |y|  -  j  <  k,  and  suppose  pg  is  expressed 

~  k  ~  ~ 

as  a  function  of  Q  for  all  s  for  which  |s|  <  J.  From  (3.5)  we  have 


JL-  £  (-i)1""-1  TT  p;1 

~  ~€lk(Z*)  Cs.e^  ~ 


M 

aw 

-  s 

i-0 


1  ^  TT 

Iv(y*)  !  |w  ■  y*!"1!  [«PA:sw*0] 

a/  R  /V  /V  w  *'  aw  as# 


-1 
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obtain 
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V  [f  M>  J  TI 


{«  :  s  y**0, | s |  <  I y| } 

M  R  MM  >v  M  0S4 


w-1  v  t 

£1  ^  ('u  tt  -s  tt 

(w€  I.  (y*)  :  |w>y*|-i}  (s  €«7.  :  |  s  J  <  |y| -1,  {s€-7.:|y|-i 

^  K  <v  /V  M  A#  K  M  M  K  <V 

$  1*1 <  W  » 

s  y*  *  0) 


(3.9) 


Equation  (3.9)  may  be  derived  from  (3.8)  by  noting  that  for  w 6  I  (y*) , 

M  K  M 

s»  ■  0  *>  s  y*  ■  0,  In  equation  (3.9),  p  Is  expressed  recursively 
as  a  function  of  Q  and  (p_  :  |s|  <  |y|},  where  1  <  |y|  <  k. 

^  8  M  ““  ^ 

M 

We  may  determine  p^  from  the  equation 

MW 

2i  ’  TT  p,  • 

•  <v  ~ 

M  K 

that  is,  as 


P,  -  Q! 


[»  EV  :  | at  <k) 


(3.10) 


*1  »  *  T 

The  estimate  p  is  obtained  by  replacing  Q  and  p  with  Q  and  p  In 

Z  ^  X>  X 

equations  (3.7),  (3.9)  and  (3.10). 
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Several  remarks  are  In  order.  First,  It  is  clear  that  £*  is  undefined 

A  A 

I 

when  nQ.  ■  N  *  0.  In  fact,  the  process  by  which  p  was  obtained 

II  am 

AM 

requires  that  p  i*  0  V  s€V.  (indeed,  this  implies  Q.  t  0  by  (3.10)). 

S  **  K  I 

AW  A  j  AW 

Finally,  even  when  p  is  well  defined  by  the  process  described  above,  it 


A  j  2  - 1 

may  well  happen  that  p  ?  (0,1]  ,  and  in  such  cases,  the  estimate  can 

am  .  V 

2-1 

clearly  be  Improved.  On  the  other  hand,  if  p  €  (0,1]  ,  then  it  is  the 


V  AM 

I  .  «  2  -1 


MLE  of  the  vector  p  ,  and  if  p  6  (0,1] 


with  sufficiently  high 


probability,  then  the  estimate  p  may  well  have  good  statistical 

AM 

properties.  In  order  to  have  a  well-defined  estimator,  let 

if  p1  €  (0, l]2  _1 


i-r. 


otherwise. 


(3.11) 


We  show  below  that  the  estimator  p  is  in  fact  asymptotically  optimal. 

AM 

*  ^  I 

(Were  jg  to  be  more  carefully  defined  when  p  exists  but  lies  outside 

(0,1]  ,  one  could  undoubtedly  obtain  an  asymptotically  equivalent 

estimator  with  better  small  sample  properties.) 

* 

Let  p  denote  the  MLE  of  p  ,  and  let  *J(p)  denote  the  information 

L  /v 

V 

matrix  of  the  MVB(2  -1)  model,  that  is, 


r 


j(p)j 


(3.12) 


for  u,  v€^.  The  computation  of  J(p)  is  straightforward  and  is  omitted. 

Theorem  3.1.  Let  Y^ , . . . ,Yn, . . .  be  a  sequence  of  iid  random  vectors 

Ic 

distributed  according  to  MVB(2  -1).  Let  p  be  the  parameter  vector  of 

aw 

the  model,  and  assume  that  p  €  (0,1)  V  y  Then 

y  ^ 
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(1)  p  p  ,  and 


(il)  (p  -  p)  — >  X~  N(0,  J'1(p))  , 

aw 

A 

that  is,  p  is  strongly  consistent  and  asymptotically  equivalent  to  the  MLE. 

/v 

If 

Proof.  The  regularity  of  the  MVB(2  -1)  model  is  easy  to  verify.  Wc 

A 

first  prove  the  consistency  of  js.  The  relationship  between  p  and  Q 
given  in  (3.7),  (3.9)  and  (3.10)  makes  clear  that  there  is  a  one-to-one 
relationship  between  p  and  Q  when  each  p  >0  and  >  0  , 

aw  a* 

Inequalities  which  hold  by  hypothesis.  Moreover,  the  relationship  is  a 
continuous  one  in  a  neighborhood  of  the  true  p  and  Q  .  Thus  p  *  f(Q), 

/v  A#  aw  a- 

A  g  g 

where  f  is  continuous.  Now  since  Q  — — ->  Q  ,  it  follows  that 

AW  AW 

f  (Q)  *•*■*>  f(Q).  Thus,  outside  of  a  null  set,  the  sequences  C  f  (Q) } 

AW  /V  Aw 

A  A 

converges  to  f(Q)  as  n  *■»  ®.  For  n  sufficiently  large,  p  «  f(Q) 

AW  AW  AW 

and  the  sequence  [p3  will  thus  converge  to  p  .  Therefore,  outside  of 

AW  AW 

-  A  a  s  A 

the  same  null  set,  .p  "•  p  and  we  may  write  p  — — ->  p.  Now  let  p 
~  ~  ~  ~L 

be  the  maximum  likelihood  estimate  of  p  .  Using  |!*j|  for  Euclidean 

AW 

distance,  we  have 


P lip  -  P  II  >  «)  <  1  -  P(P  €  (0,1)  )  -  0 

*  P  2^-1 

since  p  — >  p  €  (0,1)  .  Thus 


Since 


A  A  _ 

•Jn  (p  -  p  )  — >  0. 

AW  A/L 


Jn  (p-p)  ■  */n  (p-p)  +*/n  (p-p), 

AW  AW  awL  aw  AW  aw  Li 


we  have  that  p  and  p  are  asymptotically  equivalent,  or, 
~  ~L 


^/n(p-p)  —  >  X  ~  N(0,J_1(p)). 
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IV.  INFERENCE  FOR  SUBMODELS 


Our  treatment  of  the  general  MVB  model  is  relatively  complete.  We 
are  able  to  give  a  simple,  finite  recursive  scheme  which  specifies  a 
strongly  consistent  and  asymptotically  efficient  estimator  of  the  para¬ 
meter  vector.  Nonetheless,  the  general  model  is  unsatisfactory  in  one 
key  respect.  For  systems  with  more  than  a  handful  of  components,  the 
general  model  has  so  many  parameters  that  the  sample  size  required  for 
our  asymptotic  results  to  take  hold  could  well  be  prohibitive.  This 

difficulty  can  be  circumvented  in  situations  where  one  can  justify  the 

Ic 

use  of  a  submodel  of  MVB(2  -1)  whose  parameter  space  is  of  substantially 
lower  dimension.  In  this  section,  we  treat  the  estimation  of  the  para- 
ireter  vector  of  a  reasonable  submodel  with  a  k-dimensional  parameter  space. 

Ir 

Let  Y  be  distributed  according  to  the  MVB(2  -1)  distribution,  sub¬ 
ject  to  the  following  parametric  restrictions: 

* i-i ( \ 3  hi "  hi • 


If  1st  «  )tl  ■  J,  the  common  parameter  value  in  (4.1)  will  be  denoted  by 

p  .  We  may  rewrite  the  mass  function  of  Y  in  terms  of  the  parameters 
J  ~ 


Pj,  J*l»2,...,k,  as  follows. 


P(Y-y)  - 


?,  2  TT  77  ps 

i*|y|  ik(y>-*  j-1  ~ 

M  m  *  2 

M  mi ) 
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Using  the  convention  (  )  *  0  If  t  >  r,  the  cardinality  of  the  set 
fs€^.  :  |  s  |  =  j ,  s  w  t  0 ]  may  be  seen  to  be 

^  |v  ^  ^  ^  ^ 


k  tc-4 

<j»  *  <  j  ’ 


for  any  wfa/^  for  which  |w|  =  l.  Moreover,  the  cardinality  of  the  set 
[w€  I.  (y)  :  |w|  -£}  Is 


C  ~  ;• 

Vx-|y|y 


It  follows  that 


P(Y-y) 


,  .k.  .k-i. 

k  (,)-(,) 


which  we  rewrite  as 


P(Y-y) 


‘-•Wrk-ui'  t-i-  <$>-<k‘'-5l‘t) 


(4.2) 


The  model  whose  probability  function  is  displayed  In  (4.2)  is  a 
k-parameter  multivariate  Bernoulli  distribution  which  we  henceforth  refer 
to  as  MVB(k).  This  model  retains  some  of  the  essential  features  of  the 
general  model  we  have  studied.  In  particular,  if  Y  ~  MVB(k),  then  the 
components  of  Y  are  positively  correlated.  Moreover,  the  submodel 
preserves  the  notion  that  shocks  of  varying  gravity  may  occur  while 
making  the  simplifying  assumption  that  shocks  of  the  same  gravity  (i.e., 
simultaneously  fatal  to  a  fixed  number  of  components)  are  equiprobable. 
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One  can  easily  deduce  from  this  assumption  that  the  components  of  Y 
have  identical  marginal  distributions.  Indeed,  all  subvectors  of  Y  of 
a  fixed  dimension  will  be  identically  distributed.  The  model  should  be 
applicable  in  reliability  experiments  in  which  the  components  of  a  system 
are  inherently  similar  when  viewed  one  at  a  time  (e.g.,  the  cells  of  an 
automobile  battery),  yet  tend  to  be  positively  dependent  due  to  the 
increased  load  on  working  cells  as  each  failure  occurs. 

We  proceed  with  our  development  of  an  efficient  estimator  of  the 
parameter  vector  of  MVB(k).  Note  that  the  probability  function  in  (4.2) 
depends  on  an  observed  vector  y  only  through  !yj,  the  number  of  ones  in 
the  vector  y .  Thus,  the  likelihood  function  L  is  proportional  to 


where,  for  x-0,l,...,k. 


Q  a  P(|y|  -x)  -  L~i  P (Y  -  y) 

{y:|yj=x3 


,  k-x 

-  <k)£ 
x  t-o 


(kt*)(-D 


(k)  -  (k_x' 
v  1  j 


j 


(4.3) 


and  N  ia  the  number  of  occurrences  of  1y|  «x  in  the  sample.  If  Q  was 
x  ~  ~ 

completely  unconstrained,  the  MLE  of  the  vector  (Qq>Q^> • • • >Q^)  would  be 

/In  In  In  ),  As  in  the  previous  section,  we  will  obtain  en 
vn  O’  n  1’  n  k' 

estimator  of  the  vector  p  by  inverting  the  functional  relationship 


Q 


k 


9 


(4.3).  Since 
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we  have  (assuming  >  0) 


which  yields 


<k 


V(1/k)Vi 


(4.*) 


Now,  suppose  p^ . px_i  have  been  expressed  as  functions  of  the  Q's  for 

a  fixed  x  between  2  and  k-1.  We  then  solve  the  equation 


<k-x 


,  _  x- 1  ^  x-t  -(x”t) 

-  ( k_ )[  £  (*>(-i>t'Tr  J  + 

K  *  L  t=0  1=1  J  J 


(4.5) 


for  px<  We  note  that  (4,5)  may  be  rewritten  as 


x  -O  x-1 


►  ,x-t. 

x-t  -(  .  ) 


^-OCTTpj  j  +S(*)M)tT]'p  j  +(-dx], 

yk  •  x  Lj=i  J  t=i  c  j=i  J  J 


so  that 


<x 


x-l  x-1 

-  2  (*)(-l)CTr  P| 

t-1  J-l  3 


(4.6) 


Finally,  we  may  solve  the  equation 
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where  J(p)  is  the  information  matrix  for  the  model  MVB(k). 

The  proof  of  Theorem  4.1  is  similar  to  that  of  Theorem  3.1,  and  is 
omitted. 


In  their  paper  on  maximum  likelihood  estimation  for  the  multivariate 
exponential  distribution,  Proschan  and  Sullo  (1976)  study  in  detail  the 
submodel  in  which  the  only  shocks  which  occur  with  positive  probability 
are  shocks  fatal  to  single  components  only  or  the  universal  shock  which 


ia  fatal  to  all  components  simultaneously.  The  analogue  of  this  submodel 
In  our  context  Is  the  MVB  distribution  subject  to  the  restriction  that 

P8  -  1  ¥  s  9  1  <  |s|  <  k. 

This  is  a  k+1  parameter  submodel  for  which  we  have  obtained  the  follow¬ 
ing  results:  Direct  maximum  likelihood  succeeds  to  the  extent  that  the 
MLE  may  be  identified  explicitly  as  a  function  of  the  smallest  root  of  a 
certain  k*"^  degree  polynomial.  Moreover,  we  are  able  to  give  a  simple 
Iterative  procedure  which  converges  to  the  desired  root.  The  approach 
taken  in  our  study  of  this  submodel  differs  completely  from  that  taken 

in  this  paper.  Full  details  appear  in  Boyles  and  Samaniego  (1980). 

k 

The  two  submodels  of  MVB(2  -1)  discussed  above  serve  to  illustrate 
the  fact  that  one  may  model  the  behavior  of  positively  dependent  compo¬ 
nents  in  reliability  experiments  in  such  a  way  that  (1)  the  parameter 
space  has  a  manageable  dimension  and  (2)  efficient  estimation  of  parame¬ 
ters  la  tractable  and  computationally  feasible.  Moreover,  the  submodels 
considered  here  correspond  to  realistic  constraints,  that  is,  represent 
situations  in  which  (a)  components  seem  to  have  the  same  probabilistic 
behavior  when  viewed  one  at  a  time,  or  (b)  among  shocks  fatal  to  more 
than  one  component,  only  the  catastrophic  or  universal  shock  is  deemed 
important  or  at  all  likely. 
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V.  A  MONTE  CARLO  STUDY 


In  the  preceding  sections  we  have  derived  the  asymptotic  properties 

*  k 

of  the  estimator  p  for  the  parameter  p  of  the  MVB(2  -1)  and  MVB(k) 

***  M 

A 

models.  In  particular,  P(p  ■  MLE]  -  1  as  sample  size  -<*>.  We  now  employ 
Monte  Carlo  simulations  to  investigate  the  rate  of  this  convergence  as  a 
function  of  sample  size  and  the  true  parameter  p  . 

/ss 

Table  1  below  summarizes  the  bivariate  case  k»2  for  the  fixed 


sample  size  n*50.  For  each  parameter  vector,  P{p  ■  MLE)  is  the  observed 

A* 

A 

frequency  of  the  event  p  “  MLE,  based  on  100  simulations. 

/V 

A 

Table  1  shows  that  p  "works  well,"  i.e.,  equals  the  MLE  with  proba- 
bility  near  1,  when  shocks  are  "rare,"  i.e.,  idien  the  components  of  p 

rw 

A 

are  "close"  to  1.  In  other  situations,  it  is  evident  that  p  may  perform 

AW 

quite  poorly.  As  we  should  expect  from  symmetries  in  the  MVB  parametri- 

A 

zatlon,  the  performance  of  p  appears  to  be  invariant  with  respect  to 
permutations  of  and  p^.  One  might  also  suspect  that  p  should  work 
better  in  the  submodel  than  the  full  model,  since  the  former  has  fewer 
parameters.  However,  Table  1  does  not  support  this  conjecture.  Table  1 

A 

does  show  that  the  performance  of  p  is  sensitive  to  certain  order  re¬ 
¬ 
lations  among  the  components  of  p  .  In  particular,  p  performs  best 

when  at  least  one  of  p^  and  p^  is  greater  than  p^. 


A  A 


JT 


Table  1.  100  x  P{p  -  HLE)  (k-2,n«50). 


pl 

p2 

p10 

p01 

P11 

Pull  Model 

0.1 

0.1 

0.1 

0.1 

0.1 

7 

0.1 

0.5 

0.1 

0.1 

0.5 

17 

0.1 

0.9 

0.1 

0.1 

0.9 

37 

0.1 

0.5 

0.1 

24 

0.1 

0.5 

0.5 

65 

0.1 

0.5 

0.9 

55 

0.1 

0.9 

0.1 

34 

0.1 

0.9 

0.5 

89 

0.1 

0.9 

0.9 

72 

0.5 

0.1 

0.1 

23 

0.5 

0.1 

0.5 

73 

0.5 

0.1 

0.9 

57 

0.5 

0.1 

0.5 

0.5 

0.1 

66 

0.5 

0.5 

0.5 

0.5 

0.5 

98 

0.5 

0.9 

0.5 

0.5 

0.9 

69 

0.5 

0.9 

0.1 

87 

0.5 

0.9 

0.5 

100 

0.5 

0.9 

0.9 

95 

0.9 

0.1 

0.1 

34 

0.9 

0.1 

0.5 

92 

0.9 

0.1 

0.9 

83 

0.9 

0.5 

O.t 

88 

0.9 

0.5 

0.5 

100 

0.9 

0.5 

0.9 

93 

0.9 

0.1 

0.9 

0.9 

0.1 

98 

0.9 

0.5 

0.9 

0.9 

0.5 

100 

0.9 

0.9 

0.9 

0.9 

0.9 

100 
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A  similar  study  was  carried  out  for  the  trivariate  case  k *  3  and 
the  fixed  sample  size  n*50,  with  results  similar  to  those  of  the  pre¬ 
ceding  case.  To  save  space,  we  present  only  a  small  portion  of  this 
study  in  Table  2  below. 


Table  2. 

100  X  p{p  -  MLE] 

cv 

(k-3,n-50). 

P1 "  p100  "  p010  "  p001 

p2  *  p110  *  p101 " 

p011  p3”plll 

Full 

Model 

Sub- 

Model 

0.70 

0.70 

0.70 

82 

78 

0.70 

0.70 

0.90 

88 

68 

0.70 

0.90 

0.70 

56 

80 

0.70 

0.90 

0.90 

58 

73 

0.90 

0.70 

0.70 

98 

93 

0.90 

0.70 

0.90 

97 

77 

0.90 

0.90 

0.70 

88 

100 

0.90 

0.90 

0.90 

95 

93 

Table  2  shows  again  that  p  works  well  in  the  "rare  shocks"  region 

aw 

of  the  parameter  space,  although  a  comparison  of  Tables  1  and  2  shows 
that  the  boundaries  of  this  region  depend  on  the  dimension  of  the  model. 
Once  again  the  submodel  does  not  show  a  clear  advantage  over  the  full 
model.  The  submodel  performs  best  when  the  ordering  p1  >  p,  >  P3  obtains, 
but  the  full  model  does  not  appear  to  be  sensitive  to  order  restrictions 
in  this  case. 

A 

Tables  3  and  4  below  illustrate  the  large-sample  behavior  of  p  . 
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Table  3.  100  x  P{p*MLE]  (Full  mode  l,k*>3,  large  samples). 


Sample  Size 

P111P010P001 

P110P101P011 

Plll 

100 

200 

300  400 

500 

600 

700 

800 

900 

1000 

0.30 

0.30 

0.30 

4 

4 

8 

14 

13 

13 

17 

24 

18 

22 

0.30 

0.50 

0.70 

22 

37 

49 

55 

61 

67 

70 

69 

71 

73 

0.50 

0.30 

0.70 

23 

34 

45 

57 

68 

74 

78 

85 

88 

91 

0.50 

0.50 

0.50 

44 

82 

82 

92 

91 

94 

92 

95 

97 

99 

0.50 

0.70 

0.30 

47 

69 

80 

79 

80 

90 

88 

92 

93 

93 

0.70 

0.50 

0.30 

65 

90 

95 

99 

100 

100 

99 

100 

100 

100 

0.70 

0.70 

0.70 

94 

100 

100 

100 

100 

100 

100 

100 

100 

100 

0.90 

0.90 

0.90 

100 

100 

100 

100 

100 

100 

100 

100 

100 

100 

Table  4.  100  X  P{p*MLE}  (Submodel, k“3, large  samples). 


Sample  Size 

P1 

P2 

p3 

100 

200 

300  400 

500 

600 

700 

800 

900 

1000 

0.30 

0.30 

0.30 

1 

3 

0 

4 

12 

12 

16 

18 

18 

19 

0.30 

0.50 

0.70 

18 

21 

31 

42 

39 

53 

66 

62 

67 

70 

0.50 

0.30 

0.70 

21 

39 

38 

49 

50 

56 

64 

69 

68 

67 

0.50 

0.50 

0.50 

47 

71 

83 

89 

94 

93 

96 

93 

97 

98 

0.50 

0.70 

0.30 

69 

86 

90 

90 

91 

94 

92 

98 

100 

97 

0.70 

0.50 

0.30 

67 

89 

96 

100 

100 

100 

100 

100 

100 

100 

0.70 

0.70 

0.70 

96 

98 

100 

100 

100 

100 

100 

100 

100 

100 

0.90 

0.90 

0.90 

100 

100 

100 

100 

100 

100 

100 

100 

100 

100 

The  rate  of  convergence  of  P{p*MLE)  to  1  increases  dramatically  in 
both  models  as  we  move  toward  rare  shock  parameter  values.  The  parameter 
vectors  are  arranged  in  increasing  lexicographical  order,  and  in  both 
models  this  rate  of  convergence  is  evidently  increasing  with  respect  to 
this  ordering. 

A 

Our  final  study  involves  the  performance  of  p  for  small  and  moderate 

A 

sample  sizes.  Our  simulations  have  shown  that  here  p  will  perform  very 

IV 

poorly  except  in  the  rare  shock  case.  In  the  case  of  very  rare  shocks 

(i.e.,  all  components  of  p  greater  than  0.90)  we  found  that  the  event 
*  •** 

p  ■  MLE  obtained  in  more  than  84%  of  170,100  simulations  performed  with 
very  small  samples  (n  ■  1, 2, . . . ,7) ,  and  in  more  than  67%  of  72,900  simu¬ 
lations  performed  with  moderate  sample  sizes  (n  ■  10,20, 30) .  It  is  notable 
that  our  estimation  procedure  works  fairly  well  even  when  the  sample  size 
is  less  than  the  dimension  of  the  parameter  vectors.  The  explanation  for 

A 

this  unexpected  result  is  related  to  the  fact  that  p  appears  to  work 
better  in  very  small  samples  than  for  moderate  sample  sizes.  This  may  be 
attributed  to  the  fact  that,  for  very  rare  shocks  and  very  small  samples, 
only  a  few  data  configurations  possess  a  non-negllgible  probability  of 
occurring,  and  that  in  fact  those  turn  out  to  be  configurations  for  which 

A 

p  works.  With  moderate  sample  sizes  this  restriction  is  dropped,  and 
the  asymptotic  properties  have  not  yet  taken  hold,  hence  we  observe  a 
decline  in  performance. 


Our  simulations  have  shown  that  p  performs  very  well  in  large 
samples  whenever  the  assumption  of  moderately  rare  shocks  (e.g.,  all 
components  of  p  >  0.70)  is  appropriate.  Under  stronger  assumptions 
(e.g.,  all  components  of  p  >  0.90)  our  estimation  procedure  works  fairly 
well  at  small  and  moderate  sample  sizes.  The  simulations  also  show 
that  Che  estimation  procedure  is  more  sensitive  to  the  parameter  values 
corresponding  to  low-gravity  shocks  than  to  the  parameter  values  corres¬ 
ponding  to  high-gravity  shocks.  This  should  be  taken  into  account  when 
assessing  the  usefulness  of  our  procedure  in  a  given  situation. 

VI.  DISCUSSION  AND  CONCLUSIONS 

We  have  described  in  this  paper  a  multivariate  Bernoulli  distribu¬ 
tion  which  may  be  viewed  as  a  shock  model  in  reliability.  It  is  a  dis¬ 
crete  analogue  of  the  multivariate  exponential  distribution  of  Marshall 
and  Olkin  (1967)  but,  in  fact,  serves  as  a  discrete  version  of  any  con¬ 
tinuous  model  for  which  component  lifetimes  are  modeled  as  minima  of 
waiting  times  for  selectively  fatal  shocks.  Asymptotically  optimal 

estimators  have  been  given  in  closed  form,  both  for  the  general  model 
k 

MVB(2  -1)  and  for  a  natural  submodel  MVB(k)  with  a  parameter  space  of 
substantially  lower  dimension.  The  estimators  produced  are  only  of 
value  when  they  are  equal  to  the  MLE,  since  we  have  defined  them  quite 

A 

arbitrarily  elsewhere.  We  have  Investigated  the  rate  at  which  P(paMLE) 

AS 

tends  to  one  through  a  Monte  Carlo  study. 


The  Monte  Carlo  study  makes  the  domain  of  applicability  of  our 

A 

results  rather  clear.  We  have  found  that  P(p  =  MLE)  is  high  for  moderate 
sample  sizes  only  when  the  shocks  in  the  model  (or  its  submodel)  are  rare. 
This  has  been  found  to  be  particularly  true  for  shocks  affecting  only  one 
or  possibly  a  small  number  of  components.  The  method  of  estimation  is 
somewhat  more  robust  with  respect  to  the  rareness  of  shocks  of  greater 
gravity.  As  Beuhler  (1957)  has  pointed  out,  reliability  experiments 
often  deal  with  highly  reliable  systems  in  which  the  probability  of 
failure  of  any  given  component  is  very  low.  It  naturally  follows  that 
the  shocks  affecting  such  components  occur  only  rarely.  Thus,  the 
method  of  estimation  advanced  in  this  paper  tends  to  work  well  in  a  class 
of  problems  that  occur  frequently  in  practice. 

The  Monte  Carlo  study  has  also  revealed  that  the  rate  at  which 

A 

P(p-MLE)  -*  1  is  about  the  same  for  the  full  model  and  for  its  submodel. 

In  a  sense,  this  is  a  disappointing  result,  since  one  would  hope  to  gain 

some  advantage  when  one  dramatically  reduces  the  dimension  of  the  parameter 

space.  We  should  emphasize,  however,  that  the  rate  of  convergence  studied 

in  Section  V  does  not  discredit  the  submodel.  It  is  a  comment  only  on  the 

inversion  process  involved  in  obtaining  our  estimator.  There  is  a  second 

convergence  rate  of  interest  in  our  problem,  namely,  the  rate  at  which  the 

asymptotic  theory  of  maximum  likelihood  estimation  takes  hold.  It  is  well 

known  that  this  rate  is  affected  by  the  size  of  the  model,  and  it  is  in 

k 

this  domain  that  we  expect  the  advantage  of  MVB(k)  over  MVB(2  -1)  to  ahow 
up.  When  MVB(k)  is  appropriate  as  a  model,  we  expect  that  the  use  of 

Ir 

MVB(2  -1)  would  be  very  costly  in  terms  of  the  efficiency  of  the  MLE. 


There  are  a  number  of  Issues  that  we  postpone  for  future  Investiga¬ 
tions.  Of  particular  interest  is  the  development  of  modeling  and  infer¬ 


ence  for  larger  classes  of  MVB  distributions  with  positive  dependence. 

In  this  regard,  inference  for  the  class  of  all  MVB  distributions  with 
positively  correlated  components  may  be  feasible  using  the  general  repre¬ 
sentation  theorem  for  MVB  distributions  developed  by  Bahadur  (1961)  and 


Lazars fe Id . 
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